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Groundwater is one of the important sources of freshwater 
and accordingly, there is a need for optimizing its usage. In 
this paper. four multi-objective metaheuristic algorithms 
with new evolution strategy are introduced and compared for 


je the optimal management of groundwater namely: Multi- 
Genetic algorithms objective genetic algorithms (MOGA), multi-objective 
Memeti algoritims; memetic algorithms (MOMA), multi-objective particle 
Particle swarm: swarm optimization (MOPSO), and multi-objective shuffled 
Shulled fog leaping; frog leaping algorithm (MOSFLA). The suggested evolution 
(Compromise solution: process is based on determining a unique solution of the 
Multiple objectives Pareto solutions called the Pareto-compromise (PC) solution. 
optimization The advantages of the current development stem from: 1) 


‘The new multiple objectives evolution strategy is inspired 
from the single objective optimization, where fitness 
calculations depend on tracking the PC solution only through 
the search history: 2) a comparison among the performance 
of the four algorithms is introduced. ‘The development of 
each algorithm is briefly presented. A comparison study is 
carried out among the formulation and the results of the four 
algorithms. The developed four algorithms are tested on two 
multiple-objective optimization benchmark problems. ‘The 
four algorithms are then used to optimize two.objective 
groundwater management problem. The results prove the 
ability of the developed algorithms to accurately find the 
Pareto-optimal solutions and thus the potential application on 
real-life groundwater management problems, 
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1. Introduction 


Groundwater is considered as one of the freshwater sourc 


needed to serve human activities, 
such as drinking, industrial, domestic and irrigation. Due to irregular pumping of groundwater, 
serious environmental hazards have been occurred including groundwater levels decline and 
interference among wells 


To pump the maximum amount of groundwater without aquifer depletion and achieving 
minimum operation cost, powerful simulation-optimization models have to be applied to obtain 
the best strategy. These models are adopted to deal with both the design and operation problems 
corresponding to remediation, water supply and controlling of groundwater flow [1]. 
Optimization models can be classified according to their objectives into single objective and 
multiple objectives. In single objective optimization, a unique optimal solution is produced 
related to the maximization or minimization of a single objective [1-19]. In multiple objectives 
‘optimization, a group of non-dominated solutions, called Pareto-front (PF), are usually produced 
[20-25]. Optimization models are classified into deterministic and metaheuristic. Deterministic 
‘methods include linear, non-linear and dynamic programming [6-9,26-33]. While, metaheuristic 
algorithms are stochastic optimization techniques that mimic nature perception, learning and 
soning such as: Genetic Algorithms (GA), Memetie Algorithms (MA), Particle Swarm (PSO) 
and Shuffled Frog Leaping (SFLA) [1-5,10-25] 


Deterministic optimization methods tend to become overly complicated and require massive 
computations when dealing with complex optimization problems with numerous constraints and 
decision variables, such as real-life water related optimization problems. While, metaheuristic 
optimization algorithms are usually used due to their direct getting solutions without the need for 
gradients and initial solutions. Simulation models are coupled with optimization techniques to 
evaluate proposed objective function(s). Simulation models apply either numerical or analytical 
approaches to simulate groundwater flow and calculate the hydraulic heads of the studied 
aquifer. Numerical approaches include finite element (FEM), finite difference and boundary 
element methods while, analytical approaches use analytical element method and Fourier series 
[16,17,19} 


‘The development made in the current paper is different from previous works in many aspects. 
The use of a new evolution strategy for generating the Pareto optimal solutions in four different 
multi-objectives algorithms namely: MOGA, MOMA, MOPSO and MOSFLA. Using such 
evolution strategy resembles solving single optimization problems. Performing a comparison 
study among the four multi-objectives optimization algorithms for optimizing ground water 
management problems. While, all previous comparisons have been performed among single 
objective optimization algorithms [34,35]. 


2. New evolution strategy 


Multiple objectives optimization problems are characterized by the existence of many optimal 
solutions (non-dominated) called Pareto-Front (PF). Many different algorithms are used to solve 
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multi-objectives optimization problems. Grierson [36] proposed a methodology to determine a 
single Pareto-Compromise (PC) solution located in the PF which fairly satisfies all objectives 
This research adopts such PC solution for the fitness calculations and then driving the evolution 
process in MOGA, MOMA, MOPSO, and MOSFLA algorithms. 


The proposed new evolutionary strategy which is based on finding the PC solution is adopted 
due to their proven effectiveness in terms of the processing time and number of Pareto solutions 
obtained in the PF, as stated by Elbeltagi et al. 37]. Also, the proposed multi-objective 
evolutionary strategy is reminiscent of single-objective optimization, where its fitness 
assignment is based on tracking a single evolving solution over the search history. In addition, 
knowing the PC solution helps decision makers to choose a single solution that uniquely 
represents a mutually agreeable trade-off between all competing objectives for the problem. 


The suggested methodology for developing each algorithm to handle multiple objective 
problems depends mainly on the principle of non-dominant solutions located in more densely 
populated regions in the objective space having lower probabilities to be selected. Consequently, 
better distributions of optimal solutions can occur in the PF. The following steps illustrate the 
suggested methodology: 


1. All solutions in an evolutionary cycle are evaluated according to the values of their 
objectives to determine a set of non-dominated solutions called PF. Each non-dominated 
solution belongs to the PF can b st one objective. 

2. For each obtained PF, the PC solution is determined mathematically as given by Grierson 
[36] 

3. The distance (d,,.,) between the PC and each Pareto optimal solution (k) is calculated 
based on Eq. 1 


{all other solutions in at I 


rcs [aire 


52 Besse (o) 


ive function i 


in which, foci: the objective function i value corresponding to the PC, fi the obj 
value corresponding to solution (k) located in the PF, N: number of Pareto optimal solutions 
located in the PF, and n: number of objective functions. 


4, The fitness of each Pareto-optimal solution (&), Eq. 2, equals its distance (dy. from the 


PC solution. This step ensures the disperse of the solutions along the solution space and 
prevents solutions crowdies. 


Fitness, 


2, Beeseea 2) 


5. The relative fitness of each solution in the PF is then calculated as follows: 
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Fitness, 


>Y Fimess, 


6. To form the next generation, solutions from the current PF are randomly selected based on 
their Relative Fitness. Accordingly, Pareto optimal solutions with larger values of 
Euclidean distance have the higher probability to be selected. 


Relative Fitness 


@) 


3. Descriptions of the proposed metaheuristic algorithms 
‘A brief description of the MOGA, MOMA, MOPSO, MOSFLA algorithms are presented in the 


following subsections. 


3.1. Multiple objectives genetic algorithms 


GA simulates mechanisms of natural principle of the survival of the fittest. In GA, solutions are 
represented using chromosomes containing a group of genes, having values for the variables of 
‘an optimization problem. Each chromosome is, thes 


-valuated to determine its fitness using 
specific objective function. Then, best chromosomes are ex 
different GA operation (ic. sek 
chromosomes: 


anging their information through 


ction, crossover and mutation), as follows, to produce offspring 


= Selection: Fit chromosomes are likely to be sel 


ted to pass to next generations. 
= Crossover: Parent chromones are exchanging their gene 


information to. produce 
offspring chromosomes. 

= Mutation: It allows for random change of genetic information of individual genes. This 
pro ‘new genetic information to the current chromosomes, this may lead to 
avoid stagnation and pre-mature convergence. 


introduc 


‘The GA operators are continued for large number of generations until an optimal or near-optimal 
solution is obtained. The parameters that affect the performance of the 
cchromosom 


iA are: number of 


;, number of generations, crossover rate, mutation rate and crossover type. The 
probability of obtaining a global optimum solution is achieved by increasing the numbers of both 
population and generations but, this substantially increases processing time. 


In the proposed MOGA, for each generation, the steps from (1) to (5) mentioned in the previous 
section are carried out then, step (6) is continued to produce offspring chromosomes throu: 
crossover and mutation operators until a new population of chromosomes is formed. Also, a 
replacement strategy is carried out to replace the weakest chromosome in the current generation 
with a randomly selected one located in the PF of the previous generation. These processes are 
repeated for large number of evolutionary cycles until an optimum or near-optimal solution is 
iched. Figure | shows a flow chart of the proposed MOGA. 
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Fig. 1. Flowchart describing both the MOGA and MOMA algorithms. 
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3.2. Multiple objectives memetic algorithm 
Similar to the GA, MA chromosome’s elements are called memes, not genes. In MA, 
chromosomes are allowed to perform local search first, before being involved in the global 
evolutionary search [38]. Similar to GA, a population of chromosomes is, initially, generated 


randomly. Then, each generated solution is subjected to a local search to improve its fitne 
Afterwards, the GA operators are applied, to produ 
then performing some local search to maintain local optimality. In this paper, the local 
carried out through a pair-wise interchange (or in other words, swapping chromosome’s memes) 
as shown in Fig. 2, [38]. 


e new individuals. These new individuals are 
earch is 


chromosome AT See 
(etore local search) |_| mi x 

Chromosome A’ paleea ea 
{after local search) |_| s 


Fig. 2. Applying local search using pair-wise interchange (Elbeltagi etal. 2005), 


The maximum number of swaps is limited to 10 for each chromosome in order to maintain 
ing time. After swapping, if the chromosome’s fitness 
kept; otherwise, ignore the change. As previously mentioned, the MA param 
the parameters of the GA, in addition to the local-search mechanism 


improves, the change is 
crs are the same as 


‘The proposed MOMA is the same as the MOGA in addition to the local search previously 
mentioned. After performing each swap, a chromosome before change (chivire) is compared with 
its mate after (chy) the swaps, as Follows: 


= If the chaper dominates the chicére, then the change is kept. 

= If the chivjaw dominates the char, then the change is ignored and swaps are repeated for a 
specified number of times (10 times). 

~ Ino domination occurs, then one of the two chromosomes 
chromosome. 


randomly chosen as the new 


Figure 1 shows the flow chart of the MOMA. 


3.3. Multiple objectives particle swarm 
‘The PSO was initially developed by Kennedy and Eberhart [39]. In single objective PSO, a set of 
random particles (solutions) is initially generated to begin the evolution process. Through such 
process, three values are monitored by each particle i: its current location [X(], the best 
position it arrived [P(D], and its velocity of flying [¥(1)] where, /is the current evolution eyele. 


In each evolutionary cycle, the position of best particle (P,) is determined which is considered 
the best fitness for all particles. Consequently, the velocity of each particle is updated as given by 
Shi and Eberhart [40]: 
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U(tat)=oxv(1+c xn [P(0)- x,(D]+C, «5x12, (0)-X,(D)} 


Vu 2V,(141) 


Each particle position is updated as follows (Eq. 5) 


X,(1+1)= X,(1)+V,(0+1) (5) 


(E41): particle 7 velocity in iteration (+1); V{D: particle i velocity in iteration (1); 
i best position in iteration (1); P,(2): the global best particle position in iteration (1); 
X(D: particle i position in iteration (J); X(/+ 1): particle ‘ position in iteration (+1); @: constant 
called inertia weight; C; and C>: learning factors (C; = C2 = 2); r1 and rz: two random numbers 
(ranging from 0 to 1); Fags: the maximum value of velocity change. 


and o 


‘The main PSO parameters are: the population size; number of cycles; Vac: Cr 


In the suggested MOPSO, both Eqs. 4 and 5 are, also, adopted but, the two parameters P, and Py 
need to be determined repetitively during each cycle [41]. Initially, P; represents particle initial 
position which updated in the subsequent iterations as given by Baltar and Fontane [42] for each 
particle i 


the current 


= If the current position P; dominat. 
P, 
~ Ifthe current P; is dominated by the new position, then the new position is set the new P, 


the new position, then the new P; is set 


~  Ifno domination occurs, then the new P; is selected randomly from the two positions, 


In the MOPSO, there is no best particle position (P,) for all particles, but there are several 
equally good particles positions (PF) in each evolution cycle stored in the extemal repository. 


For each evolution eycle, the steps from (1) to (6) mentioned in the previous section are carri 
‘out to select one solution from the stored ones which is considered as P,. The flow chart of the 
MOPSO is shown in Fig. 3 


3.4. Multiple objectives shuffled frog leaping 


SELA combines the advantages of both MA and PSO algorithms. In the single objective SFLA, 
the evolution process is initialized with a number of memeplexes, each contains a number of 
frogs (individuals). Throughout each memeplex, two individuals (frogs) are determined based on 
their fitness: the best individual (Xs) and the worst individual (Xy). Also, in each cycle of 
evolution, the individual with the highest global fitness (X) is determined. Then, the individual 
with the worst fitness is improved in each cycle (not all frogs) through a process similar to the 
PSO, Accordingly, the change in worst frog position (D3) and its new position (Nirnew) are 
calculated as follows [38,43] 


Dy= Cr (Xa—Xw) CO) 
o 


Kivnen™ Xneuron + Dis 
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in which, C: the search acceleration factor; r:a random number in the range (0, 1]; and Dyae the 
maximum allowed frog’s position change 


Start 


andomly initialize all parle 
positions 


—! 


Evaluation of all objective 
functions and constraints 


Determination of the bat 
positon vector of patile 
2) 

v 
Determination of PF and 
cparation in the external 
repository 


Datearination aTTE 
compromise solution forthe 
obtained PF 


q 


Calculation of iness and 
relative fitness, Eqs. 2 and 3 
foreach particle inthe PF 


Determination of the Bat 
particle position, P, using 
oulete wheel selection 


Calculation of velocity and 
‘updating position for each 
particle, Eqs. 4 and S 


Fig. 3. Flowchart describing the MOPS 


(0 algorithm, 
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If a better individual (solution) is produced when applying Eqs. 6 and 7, the worst frog is 
replaced by the new produced one. Otherwise, Eqs. 6 and 7 is repeatedly calculated but Xy is 


exchanged with the global best frog X¢. If no improvement is occurred, a new frog (solution) is 
created randomly replacing the worst one. This process continues for a given number of 
iterations [38]. 


are: number of frogs; number of memeplexes; number of trials in each 
tions; Dyes; and C. 


In the proposed MOSFLA, there 
(Ziv) in each memeplex. Consequently, both the (Xs) and (Xy) are randomly selected from the 
corresponding solutions. Also, there is no global best frog (Xq) for all frogs, but the 
equally good frogs’ positions (PF) in each evolution cycle stored in the external repository. In 
this case, Xq is determined by applying the steps from (1) to (6) mentioned in the previous 
section in each evolution cycle. Figure 4 shows the flow diagram of the MOSFLA algorithm, 


we several solutions for (Xs) and several other solutions for 


4. Models implementation and verification 


To facilitate the implementation of the proposed four algorithms, each algorithm is coded by the 
authors using the FORTRAN language. Accordingly, all parameters related to each algorithm are 
studied. The most common parameters presented in the literature are studied here in the current 
development. Then, the proposed algorithms” performan 
objectives benchmark problems (Appendix A) where the obtained PFs are compared with the 
corresponding known true PFs of the two benchmark problems. 


is compared using two multiple 


Selection of Paramet 


values for each algorithm is an essential part of this study. As such, a 


large number of trials are performed to obtain the most suitable parameters” values that suit the 
t problems and the example application (presented later). In this step, initial param 
values are set based on relevant literature. Then, parameters’ values are changed through many 
experiments while the results are monitored. The final paramet 
model are presented in Table 1. Also, different termination criteria have been experimented with 
including the solution convergent 
experiments, the termination criterion has been set based on the number of evolutionary cycles 
{iteration-oriented) as presented in Table 1 


two 


values adopted for each 


and number of evolutionary cycles. Based on many trials and 


For the purpose of comparing the proposed algorithms performance, different metrics are used, 


namely: the generational distance (G,); the distance between the PC solution of the obtained PF 


and the true PF; the processing time; the number of Pareto-optimal solutions located on the PF 
and the mean square error between the obtained compromise solution and the nearest best 
alternative Pareto optimal solution, The values of all these metrics are calculated and presented 
in Table 2. 
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Fig. 4. Flowchart describing the MOSFLA algorithm, 
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Table 1 
Parameters’ values used in the two test problems and the example application, 
‘Algorithm Parameter names Test problem (I): ] Test problem (2) Example 
Kita function | _Kursawe function | application 
MOGA | Number of chromosomes 200 100 700 
= Number of generations: 1000 500 300 
= Crossover ratio 09 09 08 
= Mutation ratio 02 02 0.05 
= Crossover type: Uniform Uniform Uniform 
MOMA |= Number of chromosomes 200 TOO TOO 
= Number of generations: 1000 500 300 
= Crossover rate: 09 09 08 
= Mutation rate 02 On 0.05 
= Crossover type: Uniform Uniform Uniform 
MOPSO_|= Population size 100 100 700 
= Number of cycles: 500 500 300 
= The maximum velocity 035 06 5500 
«Values of learning factors 2 2 2 
Chand Ce Linear (1.4-09] | Linear [1.4.0.9] | Linear (1.3- 
= Value of inertia weight 09) 
MOSFLA [= Number of fogs 200 200 200 
= Number of memeplexes: 20 20 20 
= Number of tials before 
shuffling: 30 30 20 
= Number of iterations after 
shuffling: 100 100 300 
= The maximum change in 
frog position Dy 7 10 7000 
«The search acceleration 24 7 42 
factor C 


‘The Gy metric [44] is considered to test the performance of the developed four models for coping 
with multiple objectives problems. The Gy metric is used to measure the distance between the 
obtained and the true PFs: 


(8) 


in which, di: the Euclidean distance between the non-dominated solution i located in the obtained 
PF and the nearest one in the true PF and NV is the number of non-dominated solutions located in 
the obtained PF. A very small value for G, metric means that all the generated non-dominated 
solution nearly coincides with the true PF. 
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Table 2 
mparison among algorithms" results for the two test problems. 
‘Algorithm metric ‘Test problem (1): Kita funetion | Test problem @): Kursawe function 
Ge 0.00076 ‘00159 
Dit TIE TI96 
Processing time second OT TasT 
Moca. [Compromise sohation BIRT TH6T TR 
Test alierative Dax BID) (16.69, 4687 
SE TOOTS TOTES 
No of PF solutions TST a9 
Gr ODOR TOS 
Diss 02951 0150 
Processing time (second) 20.380 3.605 
Moma. [Compromise solution B78. 8.10) Tea) 
Test alierative GT, B09) CITE, ATI 
MSE TOOTS TOO 
No-af PF salatons 15 37 
Gr TOOOET TOORIE 
Dist O06 Corrs 
Processing time second) TSi Taz 
Mopso. [Compromise solution CATR Te, TT 
Test alierative GIS.8Iy (CITE, Ty 
TSE ToDo TORO 
No.of PF solutions 200 195 
Ge TooTos 0020 
Dist ‘O.005T 0.0301 
Processing time second) 7aTT 176 
MOSFLA, [Compromise soluion DIE RI TOIT, Ty 
Test alierative OSL SI) Te 
SE TaoNET COOOOOE 
No. of PF solutions Tor 19 


The graphical comparison between the obtained PFs, associated with each proposed model, and 
the true PF is presented in Fig. 5. It can be seen that, the obtained PFs seem regular and 
distributed throughout the solution space. The G, value is found to be 0.00076, 0.00086, 
0.00064, and 0.00102 for MOGA, MOMA, MOPSO, and MOSFLA respectively. Consequently, 
the obtained PFs are nearly similar to the true PF as per the very small values obtained of the Gu. 
Also, the minimum Gy value obtained is associated with the MOPSO algorithm. 


Figure 6, on the other hand, shows the obtained PFs for the proposed four models along with 
their corresponding true PF for the second benchmark problem. As noticed in the previous test 
problem, the obtained PFs seem distributed regularly all over the solution space. The Gi values 
are calculated as: 0.00159, 0,00051, 0.00326, and 0.00290 for MOGA, MOMA, MOPSO, and 
MOSFLA models respectively. Then, the obtained PFs are nearly similar to the true PF as per the 
very small distances of Gy metric, withthe best value associated with the MOMA. 


HA. El-Ghandour,E. Elbeltagi/ Journal of Soft Computing in Chel Engincering 2-4 2018) 01-2213 


f 


6 


40-30-20 - 00 4: 


a 
as LW 


4030-20-10 00 10 20 30 40 60 80 70 80 40-30 20-10 00 10 20 30 40 S0 G0 70 0 
Ft a 


re 
| 
rz 


Fig. 5, PFs produced by the four proposed algorithms and the true PF for the Kita function, 
The distance between the compromise solution of the obtained PF and the corresponding one of 
the true PF is calculated according to Eq, 9 and listed in Table 2 for each studied algorithm. The 
minimum distance reported is for the MOPSO with the second test problem and for the 
MOSFLA with the first test problem. 


(O} 


ras) 


Ye 


in Which, Xie and ne! the coordinates of the compromise solution corresponding to the true PF 
and Xyir and Yn: the coordinates of compromise solution corresponding to the obtained PF. 


Also, the Mean Square Error (MSE) is calculated based on the values of the eriteria for the 


PC solution and the values of the corresponding criteria f, for each of the N’ Pareto-optimal 
solutions (Eq. 10) [36]: 


MSE = (fn) xD 7 / fF 6 


on) (0) 


4 
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where, 1: number of criteria (i.e. objective functions). The solution that has the smallest MSE 
value is considered as the best alternative solution of the PC solution. 
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Fig, 6. PFs produced by the four proposed algorithms and the true 


for the Kursawe function, 


All codes are running on a personal computer of Intel (R) Core (TM) 13 CPU and having 4 GB 
Ram, The least processing time reported with the MOPSO for both test problems, The MOMA 
outperforms all other algorithms for the obtained number of solutions located on the PF for the 
second test problem, while the MOPSO outperforms other algorithms for the same metric in the 
first test problem, Generally, the results show that the MOPSO outperforms other algorithms in 
‘most of the comparison criteria as presented in Table 2. 


The results of the previous two test problems show that, the obtained PFs corresponding to the 
four models are very close to the true PFs, Also, both the rapid convergence and good results are 
obtained in the two cases. Accordingly, the proposed four models having the ability for dealing 
\with problems of multiple-objective functions. 


Example application 


In this section, the four proposed models with the new evolution strategy are applied on a 
popular groundwater management problem in order to maximize the pumping rates and 


HA, EL-Ghandour,E. Elbetagi/ Journal of Soft Computing in Civil Engincering 2-4 2018) 01-22 1S 


results is 


minimize the operation costs [24]. Also, a comparative study among the algorithms 


carried out. Both the plan of the studied aquifer and its cross-sectional elevation along with the 


required data are shown in Fig. 7. The aquifer is formed of sand and gravel with homogeneous 
and isotropic porous medium. The hydraulic conductivity (K) equals 50 m/day and it is assumed 
that the aquifer is subjected to uniform rainfall (J¥) of 0.001 m/day. The objective of this problem 


is to determine the pumping rates from the pre-specified system of the ten wells, achieving 


hydraulic constraints. The mathematical formulation of the optimization problem is as follows 
15,24} 


fr=max. [$e-rw] ay 
rnin {Eh Mla.) +a,(Q,)°(d,-h) +a,0l4, ~n)+Pti)+ Peo) (2) 
Subjected to: 
So, 2 Quy (13) 
Wy 2 Mi mins (4) 
Oran SOS Q.man | as) 
fl 4 a0 
-{, ae Nw (16) 
(yy 
(as) 


where, Q.: well # pumping rate; Nw: number of wells; PU): head violation penalty: hi, mix 
minimum head allowable at well 7 (= 0); Qs mis and Q; mac: Well { minimum and maximum, 
allowable pumping rates (equals 0 and 7000 m'/day respectively); az, a2, and as: cost coefficients 
for drilling, pumping equipment and operation of well respectively (a; = 4221 $/m, a; = 0.12 
Sim‘, a;= 0.03 S/m'); y: the penalty term; by, bs, and bs: constants indicating economies of scale 
(b1 = 0.299 and b: ); Ong: the water demand required (= 30000 m*/day); di: the depth of 
well and P(Q) : penalty term for the minimum pumping violation. 
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Fig. 7. The plan view and sectional elevation of the studied aquifer (24) 


In this application, a simulation model, based on FEM, is linked with each proposed optimization 
algorithm to predict the hydraulic head at each well location shown in the plan, Fig. 7. FEM is 
used to discretize the studied plan to number of linear rectangular elements and the two- 
dimensional flow equation, Eq. 19 is solved over each element [25] 


io) 


in which, K: the hydraulic conductivity of aquifer; ¢: the potential (= A” /2); h: the hydraulic 
head; J the uniform rainfall or uniform evaporation; Q,: the rate of injection or pumping for the 
well; (2): the Dirac delta function which equals 1 if = equal zero otherwise equals zero; and 
[Niy: number of field wells 


In th 


mulation model, Eq. 19 is written as follows: 


[4] x [pl = [F] (20) 
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in which, 4: the conductance/stiffiness matrix, (: the unknown vector of potentials and F: the 
load vector which contains the external source or sink of water and flux concentration, 


In order to make a comparison among the obtained PFs of the developed four models, in addition 
to the metrics described before, the Spacing (Sp) metric [44] is considered, Eq. 21, The Sp is used 
to measure the spread of the Pareto optimal solutions over the whole solution space. 


(a-4,) 1) 


in which, di: the minimum distance between Pareto-solution i and all other solutions located in 
the PE, Z:: the arithmetic mean of d, and NV: the number of solutions located in the PF. 
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. PEs for the example application produced by the developed four algorithms. 


A small value of Sp means that solutions located in the PF are equidistantly spaced, Figure 8 
illustrates the obtained PFs associated with each developed model. As presented in Table 3, the 
Sp metric is computed for cach algorithm as 231.6, 2842.6, 129.3, and 479.9 for MOGA, 
MOMA, MOPSO, and MOSFLA respectively. The Sp value for MOPSO model is the minimum 
then, the performance of this model is the best compared with the other models. Accordingly, the 
obtained PF from MOPSO model is more regular and distributed throughout the solution space. 
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Also, the processing time of MOPSO model (138.97 sec) is found to be the minimum compared 
with the processing time of the other models. Also, the highest number of Pareto optimal 
solutions obtained was with the MOPSO. From the visual inspection of Fig. 8, the PFs obtained 
ing MOGA and MOMA models are more uniformly distributed compared with the other two 
models. The Pareto compromise solution for each obtained PF is, also, determined as shown in 
Fig. 8. It is noticed that, the dev ompromise solutions is small. In 

mneral, the res feness of using the suggested evolution strategy in each 
model to drive the final PF for the example application. Consequently, the proposed four models, 
can be used to deal with groundwater real life applications having multiple objectives. 


Table 3 
Comparison among algorithms results forthe example application. 
Processing | Compromise TSE 
ime second) | _ sol alterna 
THAD THSOBEE, SOTTS. | TOOTS 
186891.0) 1869446) 
ERG | TAS 66 GST924, (5868.7. | OOGOOOETS 13ST 
penta 1894315) 189665.6) 
TS] TRIS SRD? 5625.6, —| OROHOSAT | 3968 
aul 189013.0) 189477.3) 
1 HOS | SORT CHOTES, CH35.1, | OOROOTTST 155 
pose: 18459) 1846634) 


In general, the MOPSO outperforms all other algorithms in most of the comparison criteria, 
which confirm with some previous studies. Also, identifying a unique PC solution helps the 


decision makers to determine a unique solution that satisfies all objectives fairly 


6. Limitations and suggested improvements 


tic models have been applied on a ground water example 
ely. Also, the models were experimented with two bench mark 
of PFs and performed well. The developed models used a unique 
solution that satisfies fairly all competing criteria of the problem "PC solution” to drive the 
fitness. As such, tracking a single point resembles optimizing a single-objective. Also, it helps 
the decision maker by presenting a unique PC solution especially in complex objective space: 


The four developed metaheu 
application and worked effecti 
problems with different shapes 


Despite their perceived benefits, the developed model still has some limitations with a number of 
possible improvements, including: 


~ For the groundwater management problem, considering other objectives in addition to the 
total pumping and total cost; 


- studying different aquifers shapes other than rectangular or near rectangular shapes.; and 


- Compare the performance of the proposed multi-objective optimization strategy with other 


strategies used in the literature 
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Conclusions 


This research applied new evolution strategy of tracking a single compromise solution to modify 
the GA, MA, PSO, and SFLA algorithms to tackle multiple objectives problems. Two multiple 
objectives standard test problems with known PFs are used to check the ability of the proposed 
models with the proposed evolution strategy to arrive at the correct PFs. The proposed four 
algorithms are then applied on a popular groundwater management example for maximizing 
‘pumping rates and minimizing operation costs. 


The suggested four models proved their ability to obtain the close true PFs for the two 
benchmark problems in addition to the groundwater application, Accordingly, they can 
effectively be used to deal with real-life groundwater management application, The performance 
of the different multi-objective algorithms is studied and compared. The MOPSO outperformed 
all other algorithms based on the used performance metrics in terms of regularity and distribution 
of PF throughout the solution space and smallest processing time when compared with other 
multi-objective optimization algorithms. In this application, a finite element-based simulation 
model, with a new idea to reduce the run time, is linked with each proposed model for evaluating 
the potential solution. 


‘Appendix A. 


A. Benchmark problem (1) 


The first benchmark problem is called the “Kita” problem with two objectives (Eq. A.1), two 
unknowns, and three constraints (Eq. 9) [44]: 


AY) 


—p fuy)=OSxty+1—p, OS yST (aay 


Subjected to: 


1 1B 

n=[cxty-—]s0, (AD) 

re[Lerr-B}e0, py a2) 
12p,+ pst Ps if PyPaP.>0 

puf2pepetps if popnps> as 
1 0 Otherwise 

where, P is the penalty term for constraints violations. In this study, different penalty functions 


are experimented with (linear, exponential, etc.) and the best results were associated with the 
proposed linear penalty function (Eq. A.3). 


‘A.2. Benchmark problem (2) 


‘The second benchmark problem is called the “Kursawe” problem with two obje 
unknowns [44]. 
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fs, ] as 


SLivew-o2 een, 


Fils,)= Soff +5sin(s,)] (as) 


‘where, the three unknowns are in the range [-5, 5] 
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